The Analytical Theory of Bulk Melting I: 
Exact Solution of the One-dimensional Atom Chain 
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We investigate theoretically the crucial role of interstitialcies that trigger the melting of a 
boundary-free crystal. Based on an interstitialcy model that resembles the J1-J2 model of frus- 
trated antiferromagnets with uniaxial anisotropy, we have calculated the exact partition function 
and correlation functions in a one-dimensional atom chain. The melting point and correlation be- 
havior of this crystal model show the applicability of Lindemann criterion and Born criterion in the 
one-dimensional case. 
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I. INTRODUCTION 

The physical picture of solid-liquid phase transition has 
emerged as a controversial issue motivated by both sci- 
ence and industry since early 20th century when diver- 
gent accounts of the instability mechanism and atom- 
scale pathway towards melting in various theoretical 
models did not seem to yield a unanimous prediction 
for the melting point of a realistic systewi^^i^. The 
widely-cited Lindemanni and BornS criteria, do not form 
an exception to this discrepancy at a glance: in the for- 
mer criterion, Lindemann proposed that melting is trig- 
gered by the avalanche of the root-mean-square atom dis- 
placement as soon as it exceeds a threshold fraction (S L ) 
of the atom spacing, where 8* L is called the critical Lin- 
demann ratio, a semi-empirical parameter initially con- 
ceived as a lattice type characteristic (but experiments 8 
suggest otherwise); in the latter criterion, Born argued 
that the vanishment of shear modulus (ReL- challenges 
this vanishment experimentally) is responsible for the in- 
ability to resist lattice destruction at the melting point. 

In this series of papers, we give a detailed presen- 
tation of how to join the Lindemann and Born crite- 
ria together. We use a model analogous to the J\-J% 
model for frustrated antiferromagnets to formulate that 
the cooperative creation and the spatial correlation of 
interstitialcies are the impetus that triggers and prop- 
agates instability in a crystal and that eventually un- 
dermines long-range order in the surface-free solid (bulk 
material). This idea is inspired by a recent molecu- 
lar dynamics simulation in which the role of intersti- 
tialcies is emphasized^, the previous understandings of 
the roles of vacancies in surface meltingil, the forma- 
tion and motion of interstitial clusters in metalsiSiiiiiii^ 
and semiconductorsiSiiLi^, early theories of interstitial- 
cies based on computer simulationsi^iSI and some ther- 
modynamic arguments2i*2£, and especially the notion of 
"virtual attraction" between defects in RefiSS. 

This paper deals with the exact solution of the 
boundary- free one-dimensional (fD) atom chain where 
two independent analytical algorithms are applied to 
evaluate the partition function. Sect. II presents a brief 



account of the methodology of our interstitialcy model; 
Sect. Ill provides an exact solution to the model based 
on the transfer matrix; Sect. IV is the summary of the 
results in Paper I; Appendix A elaborates on a solution 
to our model through Kac-Ward method, which serves as 
a valuable check to the results in Sect. III. 



II. MODEL HAMILTONIAN 

We begin our argument with the Hamiltonian of a ID 
atom chain composed of N atoms and 2N sites (See 
FIG. 1 for the site labelling): 



'tconf — 



2N 

E 

k=l 



(Jin k n k+ i + J2n k n k+2 ) 



(1) 



where the Boolean variable n k is the number of atoms 
occupying site k. n k = or I. In this paper, we as- 
sume J\ > and J2 < unless specified. In this 
Hamiltonian, the potential energy is set to be J\ for 
each nearest- neighbor (NN) atom pair, J 2 for each next- 
nearest-neighbor (NNN) atom pair and zero otherwise. 
In order to adapt to the question under investigation, we 
need two restrictions 



n k+ 2N = n k 



2N 

fc=l 



N 



(2) 



to reflect the periodical boundary condition and atom 
number conservation respectively. By using the 
transformation 2 ^!^ <j k = 2n k — 1, one could map the 
Hamiltonian above to the J1-J2 model 26 of a ID frus- 
trated antiferromagnet 



2N 

- ^2 UiOfcOfc+i + Ji^k+i) + const, 



(3) 



k=l 



where a k = <r k +2N = ±1 , and Ylk=i ° k ~ 0- When 
T = OK, sites labeled by an odd (or even) number k are 
all occupied by atoms. In the light of this, we may call 
the odd-number sites as lattice sites and call the even- 
number sites as interstitialcy sites, as shown in FIG. I, 
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FIG. 1: The way we label the lattice sites "*" and intersti- 
tialcy sites "o" in ID atom chain. 



or vice versa. The lattice sites and interstitialcy sites 
interpenetrate. 



III. THE TRANSFER MATRIX APPROACH 

A. The Partition Function for T > OK and 
Thermodynamic Properties 

To lift the constraint Y^j—i a k — 0, we introduce an 
phenomcnological external field h and consider the fol- 
lowing Hamiltonian: 



vanishes in the thermodynamic limit N — > oo (to be elab- 
orated in the next subsection), then it will be safe to say 
that Q N = Tre- w ( CT ' )/ fcBT is also the exact partition 
function corresponding to the Hamiltonian 7i can { in the 

thermodynamic limit because 5Z*=i = 0. 

In order to evaluate Qn — Tre~' H ^' 7 ' ^ kBT , we need to 
construct a 4 x 4 transfer matrix T = (t ss > ) so tha t 2 Ti 2 ? 



Qn = y^exp 



^ 2N 



4k B T 



fc=l 



X! r SlS2 T S2 s 3 • • • Tsjv-iEjv - TtT n 



(7) 



where = <J2k-i + i&2k- Such a transfer matrix must 
take the form: 



2N ^ 27V 

W (fl, ft) = ft- ^ a k + J ^ (Jl^kVk+l + J2°kOk+2) (4) 



fc=l 



fc=l 



with the condition <Jk+2N — <?k ■ If one could (i) evaluate 
the canonical partition function corresponding to 7i (a, 0) 
for any absolute temperature T > OK (to be elaborated 
in this subsection): 



Tre -W(<T,0)/fc B T 



W=±l} 



-H(a,0)/k B T 



(5) 



(fee: Boltzmann constant); and (ii) could show that the 
"ensemble average" 



1 ^ ^ 



k=l 



n{(T,o) 



n dh 



= Tr 




h=0 
H(<r,0)/fe B T 



(6) 



r 



Ti-i,— i_j Ti—ij-i+j 
2i+j,-i_j Ti+^-i+i 

e -2/3'(,7i + ,7 2 ) e -/3'Ji 

,= 2/3' J 2 



2— l+i,l— i 2— l+i,l+< 

Ti~is~i Ti-^i+i 
-M-M,i_j Ti + ^i + i 



e 2/3'(Ji-J 2 ) 
e 2/3'j 2 



(l-si)(l-x 2 ) 



1— xi 



1+XI 
1 — XI 

1 — XI 

1+XI 



1+Xl 
(1+Xl)(l-X 2 ) 
(1-Xl)(l+X 2 ) 

l+z 2 

1 — a; 2 



l+z 2 



1+xi 

1— Xl 



,2/3' J 2 

;2/3'(Jl-J 2 ) 

/ l+xi 

V i-^i 

l+z 2 
1 — x 2 

(1+Xl)(l-X 2 ) 
(1 -Xl)(l + X 2 ) 
/ l-Xl 

V i+si 



,2/3' J 2 
■/3'Ji 

■2^'(Jl + J2) 



l+x 2 
1— x 2 



/ l-xi 

V 1+xi 

/ l+xi 

V l-xi 
(l-xi)(l-x 2 ) 

(1+Xl)(l+X 2 ) 



(8) 



where j3' = (4/csT)" 1 and a;i = tanh/3'Ji, x 2 = 
tanh/3'J2. The matrix elements are obtained like the 
following example: 




= eX p{-/3'[(-l) x (1) f + (-1) x (1) J 2 + (-1) x (-1) ^ + (-1) x (1) J 2 + (-1) x (1) ^]} 

= exp(2/3'J 2 ) (9) 
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Each interior interaction (thick black edges) energy in 
the cell diagram is counted as one, while each peripheral 
interaction (thin black edges) is counted as half because 
such "bonds" are shared by two cells. 

The transfer matrix is diagonalized by T — SLS -1 , 



where L = diag [Ai, A2, A3, A4] is a matrix with four 
eigenvalues of T as its diagonal elements. Therefore, 
Q N = Af + A^ + Af + Af . Here, 



l+x%- 2x 2 x\ T (1 - x 2 ) \J{l + x 2 f 



(10) 



x 2 {x 2 x\ - 2) + x\ T (1 - x 2 ) ^{l + x 2 ) 2 x\ 



Ax 2 



{l-x\){l-xl) 



(11) 



In the thermodynamic limit, only the largest eigenvalue 
A2 survives in the final expression of log Qn — > log A2 . 
Hence we have the free energy i^onf as shown below: 

-Fconf (Ji, 'hiT) 

= Nk B T{\og [(1 - x\) (1 - xl)) - log 2 



- log[l + x 2 - 2x 2 x 1 

+ (1 - x 2 ) \j{l + x 2 f -4x 2 xl}}, 



(12) 



in the thermodynamic limit. 

From this configurational free energy F con { evaluated 
from the partition function, which is an analytical func- 
tion for all J\, J 2 when T > OK, one may draw the 
following inferences: 

(1) The occurrence of phase transition at finite temper- 
ature is ruled out because the partition function of the 
system, which reads e~ Fcan! / kBT , encounters no singulari- 
ties (as visualized in FIG. 2) in the thermodynamic limit 
when T > OK because the argument of the logarithm 
never vanishes in such cases. In other words, the system 
has only one phase for any non-zero absolute tempera- 
ture. 

(2) The configurational entropy S^cmf could be evalu- 
ated by 



conf 



1 OR 



conf 



Nk r 



Nk B dT 



= - — {(! + *§) — 

Ax 2 ksT 



[-(l + *a) 

x2 



■h 



+ (1-2:2) X 
2x 2 x\- 



ksT ksT 
(1 + x 2 y - 4x 2 x\]- 1 ' 2 } + log2 
-log[(l-x?)(l-^)] 
+ log[l + xl - 2x 2 x\ 

+ (1 - x 2 ) ^{l + x 2 f -4x 2 xf\. 



(13) 



(3) The short-range correlation functions (a k a k+ i) 
{o k a k+2 ) ((■}■ ensemble average) read: 



and 



(CfcCTfc + l) 



4k B T d 



N 



2 Fc ' 



(1 - £2)2:1 



1+2:2) - 4x 2 x\ 



(14) 



{(T k a k+2 ) = 



4k B T d 
N ~dT 2 

1 + x% - 



1 

2 Fa 



(l-x|)(l+a 2 ) 

^ (1 + X 2 ) 2 -4x 2 xj 



2x 2 



(15) 



(4) Sccmf, (<7 k <7 k+ i) and (a k a k+2 ) changes continuously 
as the temperature approaches absolute zero. 
The equation 



5 , conf(Ji, t / 2 ,T 1 )^0, asT^OK, 



(16) 



is in accordance with the third law of thermodynamics; 
the equations^ 



{a k a k+ i) 
(a k a k+2 ) 



-1, as T 
-1, as T 



OK; 
OK. 



(17) 
(18) 



agree with the intuitive picture at absolute zero: NN sites 
should contain one atom fixed at one site and one vacancy 
at the other (hence (n k n k+ i) = 0, (<T k <r k+1 ) = — 1 for all 
fc); NNN sites should be both occupied or both vacant 
(hence (n 2k n 2k+2 ) = 0, {ii 2k -in 2k+1 ) = 1, (<r k a k+2 ) = 
+1 for all k). 



B. Ensemble Expectation and Correlations via the 
Transfer Matrix 



Now we will employ the transfer matrb* 2 * 7 ' 28 to prove 
that the ensemble expectation vanishes in the ther- 



4 



modynamic limit for all k at finite temperature, which 
infers the vanishment of the right side of Eqn. @J. This 
will not only justify our procedure of lifting the "atom- 
number conservation" constraint but also show that the 
long-range order is non-existent at finite temperature. 

We find that 



= — Vs fe 

Qwfrf 



exp 



2 A 



-0 ^2 {J~L a k&k+-l + J2<3kOkW, 
fc=l 

HN 



T Sk _ lSk T, k T Sk s k+1 ) • • -Te 



(19) 



Here, A — T (d z <g> I + il <gi a z ) T, d z is the Pauli ma- 
trix and 



<8>T- 


f- il 


8> CTz 




1 + t 














1 - 


i 

























-1 - 



(20) 



Therefore, 

(S fe > = 



-Tr 



>N 



-Tr 



CT z (g)l + il(8)CT z jT ^ 

S" 1 fe® 1-Hl®^ SL* 



(21) 



As a matter of fact, 5 1 (a z ® 1 + il ® 3^)5 takes the 
shape of 



a 6 

c d 

e / 

5 ft 



(22) 



where the eight Latin letters represent the elements that 



— JM\= 2 

-J/|J 2 |=1/2 











2 3 
4A: B 77|4| 



(a) 




FIG. 2: This plots configurational free energy F con f and 
entropy £ CO nf versus temperature for various J1/J2 ratios. 
5'conf=0 at OK in all these cases. 
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(b) 

FIG. 3: (a) The inverse correlation length A -1 versus tem- 
perature for various J1/J2 ratios. Note that as T approaches 
absolute zero, so does A -1 because the correlation length is 
00 at OK. (b) The correlation function (akCk+j) versus j for 
fixed temperature and J1/J2 ratios. 
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are not necessarily zero. 



Tr 



1 + il ® a z ) SL 



- N 



= Tr 






eAf /A^ 



= 0. 



(23) 



It is now evident that in the limit that N 
which is equivalent to 



N 



(24) 



i=i 



i=i 



which is an indication of vanishing long-range order at 
finite temperature because Eqn. H24|) suggests that the 
interstitialcy sites are occupied by N/2 atoms and so are 
the lattice sites. 



From this result, we may conclude that unlike the 
short-range parameters such as (<Tk&k+i) and (akUk+2), 
the long-range order parameter £ witnesses a catastrophe 
at OK: 



1 

m 



2N 



£(-l) fe (2n fe -l) 



fc=i 



±1, 
0, 



T = OK 
T > OK 



(25) 



Therefore, there exists a solid-liquid phase transition at 
OK, a process characterized by destruction of long-range 
order and preservation of short-range order in the mean- 
time. 

To obtain the correlation between at different sites, 
we may extend the matrix method to the correlation 
function T (k, k + j) (for simplicity, j is assumed to be 
positive hereinafter): 



where 



T(k,k + j) 

((Sfc — (Sfe)) (Sfe+j — (Efc+j-))) = (E/cEfc+j) 
1 



Qn 
1 

Qn 



Tr <^ 



-Tr 



•69 



<< ; <*7 



(it) 



c/ 



dh 



da (£ 



Af 
Af 



\A 
A 3 



» 9 fe) : 



c/e *' + c?/ie A 



A 2 



(26) 



cf 



(1-3=2) 



\j {1 + x 2 f - 4x 2 x\ -xi\J~{l + x 2 ) 2 x\ - 4a; 2 



^(1 + Xa) 2 - 4X2^(1 + X 2 ) 2 - 4X2X? 



(27) 



(l-x 2 ) 



\J {l + x 2 ) 2 - 4x 2 x\ + xi \J {1 + x 2 ) 2 x\ - 4X2 



(1 + x 2 f x\ - 4aW(l + ^ - ^x 2 xj 



(28) 



and 



6 



„ = 2(log^ 



A = . log^ 



= 2 



log 



log 



l + xl- 2x 2 x\ + (1 - x 2 ) y (1 + x 2 ) 2 - 4x 2 x? 
;f - (1 - x 2 ) Xi + x 2 ) 2 x\ - 4x 2 - x 2 (2 - x 2 x^) 

1 + x\ - 2x 2 x\ + (1 - .t 2 ) \J '(1 + x 2 f - 4x 2 xf 
; 2 + (1 - x 2 ) xi ^(1 + .t 2 ) 2 2;^ - 4x 2 - x 2 (2 - x 2 xf ) 



(29) 



(30) 



With the definition of we obtain that 

(EfcEfc_|_j) = ((0'2fe-lCr2fe+2 i7 -l) - (CT2feCT2fc+2j)) 

((0'2fe-l0'2fe+2j) + (C2fc<72/c+2j-l)) 

(31) 

Therefore, (<T2k-i02k+2j-i) - {v2k02k+2j) = (transla- 
tional invariance) and 

G(0,2j + 1) + G(0,2j-1) 

= (&2k-l <J 2k+2j) + (< J 2k< J 2k+2j-l) 

i \X 2 ) i \\: 



Similarly, we explore the properties of H' k = cr 2 k-i — &2k 
(obtained by replacing the factor i in Sfc with —1) by the 
following procedures: 

(1) S^^-idz ® 1 — 1 <£> a z )S takes the shape of 



a (3 

x <5 

e 

7 1/00 



(33) 



where each Greek letter represents a non-vanishing ele- 
ment. 

(32) (2) Verify that 



= ((f 2 fe-lO'2/ c +2j-l) + (cT2k^2k+2j)) — {(<^2k-1^2k+2j) + (<T2k&2k+2j-l)) 

= 2G (0, 2j) - [G (0, 2j + 1) + G (0, 2j - 1)] 

1 s- 1 sPs- 1 SL r 



-Tr 



N 



-N-j 



Q 



(34) 



Then we combine the results above to find the recursion relations: 



G(0,2j) = 



cf 



X<t> 



dh 

i 



5rj 



(35) 



G(0,2j + 1) + G(0,2j-1) 



cf ( A3 V rf/i / A4 



(36) 



and the initial condition 



G(0,l) = (a k a k+1 ) 



(1 - x 2 )xi 



+ x 2 f - 4x2X 



(37) 



7 



which finally results in 



G(0,j) = (a k a k+j ) 

hi 1 -i- 1 *) 
+J(i+(-dO 



cf 



cf 



x<t> 



Ae 



Be 



dh 



dh 

i 



5rj 



(-If) (A'e 



B'e 



(38) 



where 



.4 



B 



A' 



B' 



cf 



(l + x 2 ) 2 (l-x 2 )-(l-x 2 ) 



xi^J(l + x 2 ) 2 x\ - 4x 2 - y(l + a^) 2 - 4x 2 x 2 



2^/ (1 + x 2 ) 2 x 2 _ 4a; 2V /(l + x 2 f - 4x 2 x 2 



dh 



-(l + ^r (l-s?) -(l-aa) 



xi^/(l + x 2 ) 2 x 2 _ 4x 2 + ^(1 + -t 2 ) 2 - 4x 2 x? 



2^ (1 + x 2 ) 2 x 2 -4x 2 \J(l + x 2 y -Ax 2 x{ 
(1 - x 2 ) 2 xi 



dh 



+ 5-q = l + 



yj(l + x 2 f x\ - 4x 2 ^J (1 + x 2 ) 2 - 4x 2 x 2 
(1 - x 2 ) 2 xi 



(l + x 2 ) 2 x 2 -4x 2 a/(1 + x 2 ) z -4x 2 x 2 



(39) 

(40) 
(41) 
(42) 



For large \j\, the e I ^ I / terms dominate the corre- 
lation function (akO'k+j), so A serves as the correlation 
length that characterizes the maximum range of informa- 
tion transfer in the ID atom chain. (This is visualized in 
FIG. 3) 

C. Discussions 

From the exact solutions of -F C onf and (£&), we find 
that the ID atom chain "melts" at OK. The catastrophe 
of long-range order and continuous change of short-range 
order compare reasonably with a real solid-liquid phase 
transition in three dimensions. This melting point in ID 
atom chain is in perfect agreement with the Lindemann 
criterion in the ID case, recalling the well-established re- 
sult that the mean displacement of a ID elastic atom 
chain is oo for any positive temperature 32 . Although in 
the model Hamiltonian Eqn. (JTJ, we do not explicitly take 
vibrational energy into account, the melting point in our 



model still coincides with the prediction based on vibra- 
tional instability for two reasons: (1) the pair interaction 
mode in our model (a repulsion Jj for NN, an attraction 
J 2 for NNN, say) still forms a reasonable caricature of 
the potential well in a solid. (2) the "hopping" between 
lattice sites and interstitialcy sites, which is the mode 
of atom movement in our model, mimics the quantized 
motion of atoms at low temperature. 

One may also check Born's scenario in our ID exact 
solution. In principle, there is no such a counterpart of 
the three-dimensional shear moduli in the ID case. Nev- 
ertheless, for fixed J 2 (bonding energy), greater J\ infers 
greater energy gap between the NN contact and NNN 
contact modes. In other words, when the bonding energy 
J 2 and temperature T is fixed, increasing J\ will add to 
the difficulty of creating interstitialcies, thereby increas- 
ing the "rigidity" of the ID atom chain. In the light of 
this, we may define a dimensionless "rigidity parameter" 
as ( Ji + | J 2 |) / | J 2 \. It can be verified analytically that 
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dS con f ( Jl, J 2 , T) 

(l-z?) 2 (x 2 - 1) (x 2 + if 



(1 - X?) / aS'conf (X U X 2 ) 



4k B T V 



9xi 



4fc B T /7 I : 

J[(l + x 2 ) -±x\x 2 

< 0, (0 < -x 2 ,:ri < 1) 



(tanh 1 x\ — 2xi tanh 1 x 2 ) 



(43) 



and 



(d\{Jx,J 2 ,T) 



(d\( Xl ,X2) 



A 2 2(1 -3:2) 



4fc B T 



:? (1 + x 2 f - 4x2 v /^(l + ^ 2 ) 2 -4x 2 _ 



> 0.(0 < -x 2 ,xi < 1) 



(44) 



Physically speaking, these two inequalities are reasonable 
because the increase of rigidity will help to restore some 
short range order to the system, thereby lowering the en- 
tropy and increasing the correlation length. (See FIG. 2 
and FIG. 3(a) for specific examples.) 

The exact solution above provides insight for the 3D 
case in that correlation between atom occupancies stands 
to be the common means to propagate instability in sys- 
tems, no matter what dimensionality. In a proceeding 
paper, we will present a detailed analysis of how such cor- 
relation undermines the long-range order in a fee crystal. 
Some numerical corollaries of this ID model may also 
shed light on the melting point formula to be established 
in the 3D casei^ One heuristic approach to "guess" the 
melting point formula in 3D is to find the "critical" tem- 
perature T c in the ID model where the heat capacity has 
a maximum, and speculate that such a T c will correspond 
to a heat capacity discontinuity in the 3D case. A lengthy 
expression for the heat capacity 



c = T- 



,dS 



dT 



(45) 



in our ID model could be obtained analytically. FIG. 4 
presents a "phase diagram" for the sign of dc/dT based 
on the analytical expression, in which the phase bound- 
ary (where dc/dT = 0) could be approximated by a sim- 
ple formula^: 



tanh ■ 



U 



tanh • 



U2 



0.. 



(46) 



When |Ji| 



Ak B T c 4k B T c 
I J 2 1 , the formula above could be replaced 



by 



U| + |J 2 | -3.52k B T c 



(47) 



without loss of precision. This suggests that in the sce- 
nario of 3D melting at a temperature corresponding to 



-0. 2 



dc/dT<0 



c 

CD -0. 6 



-0. 8 



dc/dT>0 



0. 2 



0. 1 



0. 6 



0. x 



tanhflJJ/4/CeT) 



FIG. 4: A "phase diagram" for dc/dT (c: heat capacity; T: 
absolute temperature) in which the phase boundary is ap- 
proximately a straight line. 



such a T Cl the thermal energy ~ 3k B T c is roughly tanta- 
mount to the energy barrier of interstitialcy creation, or 
in other words, the "rigidity" is combatted by thermal 
motion. As long as such a positive T c is existent, we hold 
the following to be true: 



tanh ■ 



Ua 



4fc B T c 



<0. 



,T 



> 



Ua 



(48) 
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IV. SUMMARY 



We find an exact solution to the ID problem, in which 
the atom chain "melts" at OK, which is consistent with 
the Lindemann criterion in the ID case. The behavior 
of the atom chain also suggests that softer lattice has 
greater disorder, because for fixed temperature T and 
"bonding energy" | J2I, both the configurational entropy 
(FIG. 3(a)) and the inverse correlation length (FIG. 3(b)) 
increase monotonically with the decreasing of the dimen- 
sionless "rigidity parameter" (Ji + \J 2 \) / \ J%\ - manifest- 
ing the reasonability of Born's criterion. 
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term of the polynomial which contains all the variables 
(Tfe in even powers gives a contribution to the sum which 
is proportional to the number of configurations, 2 2N . 

(2) We notice that each term of the polynomial can 
be uniquely related to a set of lines or "bonds" join- 
ing various pairs of adjoining lattice points as shown in 
FIG. (lA4l) . 



k-2 k k+2 k+A k+6 k+8 k+W 




k-1 k+l k+3 k+5 k+7 k+9 (A4) 

For instance, the diagrams in FIG. (|A5|) correspond to 
the terms 

(a) x 1 x 2 o k +\o\o k + 2 , 

(b) xfxlal_2(T 2 k _^al +1 al +2 ^ 

(c) — x 1 x 2 a k _ 2 a k _ 1 <r k+1 a k+2 a k a k+ll a k+3 a k+5 a k+6 . 



APPENDIX A: AN GRAPH THEORETIC 
APPROACH IN ID ATOM CHAIN 

An Exact Graph Theoretic Approach for ID 
Chain Composed of Finite Atoms 



The main idea of our graph theoretic solution is par- 
allel to that of Kac, WardM, and Vdovichenkc 34 , as re- 
ferred to in some textbooks, (e.g. See Refi^LS) However, 
the Feynman rules involved in our calculation are non- 
trivial and should be carefully developed. Furthermore, 
the evaluation of the thermodynamic limit requires con- 
siderable techniques. For the completeness of this Ap- 
pendix, we will elaborate on the solution process in three 
steps. 

(1) We use the identity 



exp(r6») = cosh6»(l + rtanh6») , r = ±1 



(Al) 



and substitutions x\ — tanh( Ji/4/csT), x 2 — 
tanh( J 2 /AkgT) to transform the partition function into 



= 11- xi 



where 



, -N 



2N 



1-xi 



, -N 



S (xi,x 2 ) 



(A2) 



S(x ll x 2 )= 2J H (! _ xia k <Jk+i) (1 - x 2 a k a k+2 ) 

{a k =±l}k=l 

(A3) 

The summand in Eqn. (|A3I) is a polynomial in the vari- 
ables x±,x 2 and a k . Each a k can appear in the polyno- 
mial in powers ranging from zero to four. After sum- 
mation over all the {a k = ±1} states, the terms contain- 
ing odd powers of a k would vanish. Therefore, a non- 
zero contribution comes only from terms containing a k 
in powers of 0,2 or 4. Since <r k = a\ = at = 1, each 





k+2 



k+A 



k+6 




k+l 



k+?, 



k+5 



(A5) 



Each level line (alias "edge") in the diagram is assigned 
a factor (—x 2 ) while each oblique line a factor (— x-\). 
Each end of each line is associated with a factor a k . 
The fact that a non-zero contribution to the partition 
function comes only from terms in the polynomial which 
contain all the a k in even powers signifies geometrically 
that either 2 or 4 bonds must end at each point in 
the diagram^ Hence the summation is taken only over 
closed diagrams, which may be self-intersecting (as at the 
point numbered k in FIG. (|A5I) (b)). 

In the light of this, the sum S {x\, x 2 ) may be expressed 
in the form 



S (xi,x 2 ) 



2 2N Y,(-^ r 9r(^) (A6) 



where g r is a polynomial with argument x\jx 2 . This 
polynomial represents the contributions from closed dia- 
grams formed from a number r of bonds, with the coef- 
ficient of {x\/x 2 ) q being equal to the number of distinct 
r-bond long closed diagrams with an (even) number q of 
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oblique lines, and each multiple diagram (e.g. FIG. (|A5|) into one over all possible loops, then we calculate it by a 
(c)) being counted once. "random walk" method and arrive at the following exact 

(3) We convert the summation over closed diagrams formula for < — x%, x\ < 1: 



The 



S (X!,X 2 ) 




2 2N 




2 2N -l 
n=0 


(2N-1 r 

I fe=0 L 


C 27V 

n[((H 


-4) (i 



l — Xl ( 1)( L 2 fe mod 2ivj+[ 2 (fc + l) ' 

\ 1/4 



1 + xl) (1 - x\) - 2x 2 (1 - xj) f (p)Y + 2x\ (1 " X2T (1 + / (p)) 



(A7) 



where 



/ (p) = cos 



(2p + l - (iV mod 2))tt 
TV 



[xj = the least integer that is greater than or equal to x. 

Remark: Every factor in the product above is non- 
negative in that it reaches the minimum value of 
(l + a; 2 ) 4 (l-x2) 2 >0at/(p)=-l. 

We will regard each diagram as consisting of one or 
more closed loops. For non-self-intersecting diagrams 
this is obvious; for example, the diagram in FIG. (|A5|) 
(c) consists of two loops. For self-intersecting diagrams, 
however, the resolution into loops is not unique: a given 
diagram may consist of different numbers of loops for 
different ways of construction. This is illustrated by 
FIG. (|A8|) . which shows three ways of representing the di- 
agram in FIG. (|A5(I (b) as one or two non-self- intersecting 
loops or as one self-intersecting loop. (When referring to 
the number of intersections, we should be heedful that ev- 
ery turning point should be "polished" before being taken 
into computation.) Any intersection may similarly be 
traversed in three ways on more complicated diagrams. 



single contribution to the sum, as they should. The new 
sum will also include diagrams with "repeated bonds" , 
of which the simplest example is shown in the far left of 
FIG. lAll) : 






(A9) 

These diagrams are not permissible, since some points 
have an odd number of bonds meeting at them, namely 
three, but in fact they cancel from the sum, as they 
should: when the loops corresponding to such a dia- 
gram are constructed each bond in common can be tra- 
versed in two ways, without intersection (as in the mid- 
dle of FIG. ilA9(l ) and with self-intersections (far right of 
FIG. HA9|l ): the resulting set of loops appear in the sum 
with opposite signs, and so cancel. We can also avoid 
the need to take into account explicitly the number of 
intersections by using the following Feynman rules: 

Feynman rules for vertices 




(-) 



(") 



(+) 



(A8) 



It is evident that the sum l|A6(l can be extended to 
all possible sets of loops if, in computing the contribu- 
tion of diagrams to g r (xi/x2), each diagram is taken 
with the sign (—1)™, where n is the total number of self- 
intersections in the loops of a given set, since when this 
is done all the extra terms in the sum necessarily cancel. 
For example, the three diagrams in FIG. (| A8|> have signs 
+, — , + respectively, so that two of them cancel, leaving a 




(AlO) 
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Feynman rules for edges 



(All) 



It is a direct corollary from the rules above that every 
four-vertex contributes a factor 1 in general. It can be 
verified that when these rules are applied to diagrams 
in FIG. l|A8p . the product of all contributions from the 
edges and vertices results in the sign as (— l) s where 
n is the sum of number of self-intersections from a set 
of s loops^ Here is another more sophisticated example 
where contributions from the six oblique edges are not 
explicitly marked: 



1 



-1 



1 



-1 



(A12) 

It could be verified that the contribution of a directed 
graph does not change if all the arrows are reversed. An- 
other thing that should be considered is that loops that 
"wrap up the circle" should bear a correct sign. When 
N is odd, this is automatically sufficed when the Feyn- 
man rules above are applied. For instance, the loop that 
connects all the consecutive sites together (and connects 
"2iV" to "1") in a circle has the sign " — ". However, 
when N is even, we cannot achieve a correct counting 
unless we impose the following additional Feynman rules 
to the directed edges: 

Additional Feynman rules for edges when N is even 



(A13) 

where e — tt/2N. 

We use these ad hoc Feynman rules instead of the "ge- 
ometric result about the total angle of rotation of the 
tangent vector" (a variant of Umlaufsatz in differential 
geometry) as did in Ref. (— ), because the latter method 
would result in different recurrence relations (see the ar- 
gument below) for odd number sites and even number 
sites, which would plague the diagonalization procedure. 

Let f r {xi/x 2 ) denote the sum over single loops of 
length r (i.e. consisting of r bonds), each loop carry- 
ing a factor ±1 as determined by the Feynman rules, and 
the power of X\jx 2 denoting the number of oblique edges 
in that loop. Then the sum over all pairs of loops with 
total number of bonds r is 



1 

2! 



f-^ \x 2 J \x 2 

n+r2— r x / x 



(A14) 



the factor 1/2! takes into account the fact that the same 
pair of loops is obtained when the suffixes t\ and r 2 are 



interchanged, and similarly for loops of three or more 
loops. Thus the sum becomes 

S (X!,X 2 ) 



2 2N 



oo 1 oc 

= £(-d ^ E (^) ri 



H Hs 



s=0 



ri,r 2> — =1 



frA- )-M 

x 2 \x 2 



(A15) 



Since S (xi, x 2 ) includes sets of loops with every total 
length r*i + r 2 + • • • , the numbers 7"i , r 2 , • ■ • in the inner 
sum take independently all values from 1 to oo. Hence 



E, \riH \-r e f ( x l 
(-X 2 ) f n 

n,r 2 ,-=l v z 



frA- 



and S (xi,x 2 ) becomes 
S (xi,x 2 



(A16) 



)2iV 



eX P I - E ( _X2 ) r f r 




(A17) 



It is now convenient to assign to each lattice point 
(marked by an asterisk in the figure below) the four pos- 
sible directions from it and to number them by a quantity 
v = 1, 2, 3, 4, say, as follows: 




(A18) 

We define an auxiliary quantity W r (k 7 v) the sum over 
all possible paths of length r from some given point fco, 
to a point k, v (each bond and vertex having as usual the 
factor contributed by Feynman rules); the final step to 
the point k, v must not be from a site marked by the 
left-right mirror of v. (e.g. 1 ^ 4, 2 <-> 3) With this 
definition, W r (ko,uo) is the sum over all loops leaving 
the point ko in the direction vq and returning to that 
point. It is evident that 



fr 



lj2W r (ko,u ): 



(A19) 



both sides contain the sum over all single loops, but 
^2 W r contains each loop 2r times, since it can be tra- 
versed in two opposite direction and can be assigned to 
each of r starting points on it. 

From the definition of W r (k, v) we have the recur- 
rence relations (sometimes referred to as Chapman- 
Kolmogorov equations) when N is an even number: 
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W r+1 (k, 1) = -e 2ie W r (k - 2, 1) + i^e ie W r (k-1,2 
W r+1 (k, 2) = e 2l£ W r (k - 2, 1) - i^e i£ W r (k -1,2)- 

+ - if-e~ ie W r (k 
1,2) 



W r+1 (k,2) = e 2l£ W/ r (fc 
W r+1 (fc, 3) = -ie 2l£ W r (k - 2, 1 
W r +i (fc, 4) ^0-fie l£ ^ r (fc 



i|i 2 e -^iy r (fc + l,3) 



f^e- l£ iy r (fc + l,3) + 
2ie W r (fc + 2,4) 
-°- 2l£ W r (fc + 2,4) 



+ ze 
1,3) 



(A20) 



e^ 2l£ VF r (fc + 2,4) 



The method of constructing these relations is evident: 
for example, the point k, 1 can be reached by taking the 
last (r + l)th step from the directions marked by 1,2,3, 
but not from the right (direction 4); the coefficients of 
W r arise from the Feynman-rule contributions of the last 
bond and vertex encountered. We can verify that these 
coefficients are the same for odd number k and even num- 
ber k. Let A denote the matrix of the coefficients in Eqs. 
1201) (with all k), written in the form 



Wr+l (k, V) 



E 

k' 



K(kv\k'v')W r {k',v'). 



The method of constructing these equations enables us to 
associate with this matrix an intuitive picture of a point 
moving step by step through the lattice with a "transi- 
tional probability" per step from one point to another 
which is equal to the corresponding element of the ma- 
trix A. The point traverses only one bond per step. It 
is evident that the "probability" of traversing a length r 
will be given by the matrix A' r . In particular the diag- 
onal elements of this matrix give the "probability" that 
the point will return to its original position after travers- 
ing a loop of length r, i.e. they are equal to W r (k Q ,u Q ). 
Hence 



TrA r = Y, W r (k ,v ), 



ka,v 



where the Xi are the eigenvalues of the matrix. Substi- 
tuting this expression in Eqn. (|A17|) and interchanging 
the order of summation over i and r, we obtain 



S (xi,x 2 ) 

2 2N 



exp | i log (1 + x 2 Ai)| 



(A21) 



The matrix A is easily diagonalized with respect to the 
suffix fc by using a invertible transformation: 



W r {p,v) 



2N 



? ipkE Wr{k,v) 



(A22) 



where e = w/2N. Taking "Fourier components" on both 
sides of Eqs. (|A20|) . we find that each equation contains 
only W r (p, v ) with the same p, so that the matrix A is 
diagonal with respect to p. For a given its elements 
are 



A {pv\pv') = 



le 





-(2+4p)ie 



_ e (2+ip)ie j£i e (l+2p)ie + £i & -(l+2p)ie 

e (2+ip)ie _^|i e (l+2p)ie 2 q 

_j e (2+4p)ie 2 q _^£i e -(l+2p)i£ -{2+ip)ie 

_£l e (l+2p)»e j£L e -(l+2p)ie „-(2+4p)i£ 

x 2 X2 



(A23) 



For a given p, a simple calculation shows that 



Y[ (1 + x 2 XJ = det ((W + x 2 Ku') 



= (l + x 2 2 )(l-x 2 1 )-2x 2 (l-x 2 1 )co S {2p+1)7T - 



N 



2ix± (1 — x%) cos 



{2p+ 1)tt 
2N 



(A24) 



where A w is the eigenvalue of block A (pv\pv'). Hence, 



2iV 



A (pv\pv) = A ipu\pv') 

p=i 
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infers that 



S (x 1 ,x 2 ) 



127V 



' 2N r 

n 

kP=1 



(l + xj) (l - xf) - 2x 2 (l - x\) COS 



(2p+1)tt a (2p+l)7r" 

— llX\ (l — X2) cos ■ 



1/2 



AT 



2N 



.(A25) 



For the odd AT case, simply replace 2p+ 1 with 2p and this 
would yield the correct result. The product in Eqn. (|A7|) 
is actually obtained by noticing that f r (—xi/x2) = 
f r (+xi/x2), and by multiplying the Eqn. (|A25ll by its 
complex conjugate. 



integral. It can be seen that both odd N and even A" 
lead to the same limit in the following integral, as they 
should: 



The Partition Function in the Thermodynamic 
Limit 



In the thermodynamic limit, when N tends to infinity, 
the summation involved in the logarithm of the partition 
function can be legitimately replaced by a corresponding 



^2N l ° gQN 



"2 [Ml 



4) 



where 



log (l 



!og2 + — I(xi,x 2 ) 



271 



I(x u x 2 ) = J log{((l + ^) (l-xf) -2a; 2 (l-2;2)cos20) 2 + 2a;2(l-2;2) 4 (l + cos2(/ ) )}d0 (A26) 

The integral involved in the partition could be evaluated with some complex analysis techniques as follows: 
Lemma: Assume that a > 0, log 1 = and that t, s are real numbers, we have 

r2-rr 

h(a,t,s) = / \og[a(l + t 2 +s 2 + 2t{l + s) cos m(j) + 2s cos 2m<m I d(f> 
Jo 

= 27rloga 



where 



Proof: 



m = 1, 2, ■ • • , min |l + tz + sz 2 \ > 0. 

MO 1 1 



1\ (a,t, s) — 27rloga = 2Re <x> log (l + tz 

J\z\ = l 



dz 



= 2Re 27rires 



,0 =0 



iz 



(A27) 



(A28) 



Corollary: 



/ log (a (l + i 2 + 2icosm0)) d0 = 27rloga,a > 0, -1 < t < 1, m = 1, 2, • 



(A29) 



With these two lemmas, we are able to set out to evaluate the integral: 



h (..",../■_> i / log{((l + ^)(l-^)-2a; 2 (l-^)cos2< ? !>) 2 + 22:2(l-a; 2 ) 4 (l + cos20)}' 

log{(l + x 2 ) 2 (1 - x 2 ) 2 + 2x\ (1 - x 2 f + 2x\ (1 - x 2 ) 2 
+ (-4x2 (1 + ^2) {l-x 2 ) 2 + 2x\ (1 -2>>) 4 ) cos2(j) + 2xl (l - x\) 2 cos40}d0 



(A30) 
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and reach the conclusion that 
Theorem: 



h (xi,x 2 ) 



'2 71 log <{ i + \<c\ + 2x\x\ + x\x\ - 2x 2 x 2 - 2x 2 x 2 



+ i (1 - x 2 ) (x\ - 2xlx 2 + 1) \J {1 + x 2 ) 2 - 4x\x 2 



Proof: With 



and 



< -x 2 < 1,0 < xi < 1 



A = (l + x 2 ) 2 (l-x 2 ) 2 + 2x 2 (l-x 2 ) 4 + 2x^(1 
B = -4x 2 (l + x 2 ) (l-x 2 ) 2 + 2x 2 (l-x 2 ) 4 
C = 2x 2 (l-x 2 ) 2 



(A31) 



(A32) 



we obtain 



a (1 + t 2 + s 2 ) = A > 
2ai(l + s) = B > 
2as = C > 0, 



a(l + s + ty = A + B + C 
a(l + s-t) 2 = A-B + C 
s C 



(l + s)i S 



(A33) and further that 

(A34) 

(A35) 



(A36) 
(A37) 

(A38) 



(1 + s) 2 _ B y/A + B + C + y/A-B + C 
s ' ~ C yj A + B + C - yJA-B + C 



(A39) 



C_ 



4 (y/A -B + C - yJA + B + C) 



1 1 

2 + 2 
> 



\M - S + C (S + 2C) + yj A + B + C (B - 2C) + V2y B 2 (a - 3C + \M - B + Cy/A + B + c) 
+ ^ + 2x 2 x 2 + xrjx 4 - 2x 2 x 2 - 2x\x\ + ]- (1 - x 2 ) (x| - 2x 2 x 2 + l) ^(1 + x 2 ) 2 + 4x 2 - x 2 



r 



(A40) 



From this, it is evident that s (xi, x 2 ) = C/2a (x\,x 2 ) 
is continuous with respect to (x\,x 2 ) and s (xi,x 2 ) > 
for < -x 2 < 1. 

We are able to affirm that I 2 (x\,x 2 ) — 2tt log a (xi, x 2 ) 
(0 < — x 2 < 1,0 < x\ < 1) if we can further that 
min^|^i |l + t (xi, x 2 ) z + s (x\,x 2 ) z 2 \ > is true for all 
{x\,x 2 ) such that < — x 2 < 1,0 < x-y < 1. To ver- 
ify this, we notice that for < — x 2 < 1, a(0,x 2 ) = 
l,t(0,x 2 ) = — 2x 2 , s (0, x 2 ) = 2x 2 , so it is easy to find 



mini. 

2 



1^! |l + t(0,x 2 )z + s(0,x 2 )z 2 | = min^i^j 

> 0. In other words, the roots of the equa- 
tion l + tz + sz 2 lie outside the unit disk (\z\ < 1) when 
— x 2 = 1/2, xi = 0. For an arbitrary point (xj,x 2 ) that 



that 

|1 - x 2 z| 



satisfies < 



< 1,0 < xT < 1, we may connect it to 



the point (0, 1/2) with a line segment. It is obvious that 
the two roots of 1 + t {x\, x 2 ) z + s (xi, x 2 ) z 2 vary con- 
tinuously as (xi,x 2 ) moves along this line segment. The 
orbits of the two roots must be two continuous paths 
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in the complex plane in the process mentioned above, 
and each of them begins with a point outside the unit 
disk and cannot end up with a point inside the unit disk 
unless it hits the unit circle (\z\ = 1) some time. How- 
ever, neither root can hit the unit circle because it is 
true that mino^e<27r |l + 1 (xi, X2) e l6 + s (xx, X2) e 2ie | = 
(1 + x 2 ) 2 (l — xf) / y/a {x\,x 2 ) > 0. Therefore, both 
roots of the equation 1 + 1 (ar*, x\) z + s (a;*, x%) z 2 should 
lie outside the unit disk. This completes the proof of 
h (xi,x 2 ) = 2nloga(xi,x 2 ) (0 < -x% < 1,0 < x\ < 1), 
and the generalization to x 2 — or xi — case is a trivial 
calculation. 



It can be verified by simple algebra that Eqn. (|A31|) 
gives the free energy expression that is exactly equivalent 
to Eqn. I|12|) . This circuitous approach provides an inde- 
pendent method to evaluate the partition function in the 
ID atom chain. 
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